function [alp,lan,mug,ome,fitted]=fitcumnomalizedcreationrate(xdata,ydata)
% 
% Usage : [alp,lan,mug,ome,fitted]=fitcumnomalizedcreationrate(Xobs,Yobs)
% 
% Matlab function from  Taisne and Tait "The effect of solidication on a propagating dike", JOURNAL OF GEOPHYSICAL RESEARCH 2010.
% 
% Input : xdata - normalized values of "creation" rates. "Creation" being a parameter describing intrusion dynamic
%         ydata - normalized cumlate numbers of each "creation" rates. See the "hist" command
%         
% Output : alp - can be seen as a proportion of values between the intrusion being halted and undergoing propagation
%          lan - free fitting parameter
%          mug - corresponds to the normalized mean creation rate during propagation events but also to the mean value on histogram
%          ome - standart deviation of histogram values
%          fitted - fitted values of Yobs
%
% EXAMPLE :        
% typical inputs are generated by the commands below ("creations" is the ocurence time of the parameter discribing intrusion dynamic)
%         load('clusterEQ2008dec.mat'); % for example
%         Cumulated = (1:length(creations))/length(creations);
%         Normalized = (creations-min(creations))'/(max(creations)-min(creations));
%         figure('position',[2 490 1392 294]);
%         subplot(1,4,1);plot(Normalized,Cumulated);axis tight;title('describing parameter');drawnow ; pause
%         diffs = diff(1:length(creations))./diff(creations);
%
% %The following is sometime necessary to get ride of the "clipped" values : 
%         diffs = diffs(diffs<median(diffs)+std(diffs)/2) ;
%         maxi = max(diffs);
%         [n,diffbins] = hist(diffs,[0:maxi/70:maxi]); 
%         ax(1)=subplot(1,4,2);hold on;bar(diffbins,n);axis tight;title('Histo of the creation rates');drawnow; pause
%         normalized_N = cumsum(n./sum(n)); 
%         nomalized_Diffs = diffbins./max(diffbins) ; 
%         ax(2)=subplot(1,4,3);hold on;plot(nomalized_Diffs,normalized_N,'-g','linewidth',3,'parent',ax(2));title('Normalized and cumulate values of previous histogram'); drawnow; pause
%
% %Then the function fitcumnomalizedcreationrate.m can be called :
%         [alp,lan,mug,ome,fitted]=fitcumnomalizedcreationrate(nomalized_Diffs,normalized_N);
%         axes(ax(2));hold on;plot(nomalized_Diffs,fitted,'--r','linewidth',2,'parent',ax(2)); legend('data','fit');drawnow; pause;
%
% %Finally dimensionless flux and temperature domains could be plot by the command below, it also output the intersection values :
%         ax(3)=subplot(1,4,4); dimlessTempNFlux(alp,mug); % Read Taisne and Tait "The effect of solidication on a propagating dike" JOURNAL OF GEOPHYSICAL RESEARCH 2010 
% 
% Fred.Massin@utah.edu, 2010


alp= linspace((0),(1),2) ; % [0 1]
mug= linspace((0),(1),2) ; % [0 1]
lan= linspace((0),(1000),2) ; % >0
ome= linspace((0.00001),(1000.00001),2) ; % >0

cnt = 0 ; 
pts = 0 ;
rms = 0 ; 
params = [0 0 0 0] ; 
%x = zeros(length(xdata),length(alp)*length(mug)*length(lan)*length(ome));
x = zeros(length(xdata),1);
y=x;
chi2m=1000;
lbo = [0;0;0;0.00001];
ubo = [1;1000;1;1000.00001];
options.Display='off';
options.Diagnostics='off';
options.FunValCheck='off';

for a=alp
    for m=mug
        for l=lan
            for o=ome                
                [par,resnorm] = lsqcurvefit(@myfun,[a l m o],xdata,ydata,lbo,ubo);                
                if (resnorm<chi2m)
                    chi2m=resnorm ;
                    parmen = par ;
                    cnt=cnt+1;
                    params(cnt,1:4) = par([1 3 2 4]) ; %[a m l o] ;                    
                    x(1:length(xdata),cnt)=xdata; 
                    [y(1:length(xdata),cnt),leg] = cumnomalizedcreationrate(xdata,par(1),par(2),par(3),par(4)) ;
                    leg = ['rms=' num2str(resnorm) ' ' leg] ;
                    legen(cnt,1:length(leg))=leg;
                end
            end
        end
    end
end
for iteration = 1:5    
    [par,resnorm] = lsqcurvefit(@myfun,parmen,xdata,ydata,lbo,ubo);        
    if (resnorm<chi2m)
        parmen = par ;
        chi2m=resnorm ;
        cnt=cnt+1;
        params(cnt,1:4) = par([1 3 2 4]) ; %[a m l o] ;        
        x(1:length(xdata),cnt)=xdata;
        [y(1:length(xdata),cnt),leg] = cumnomalizedcreationrate(xdata,par(1),par(2),par(3),par(4)) ;
        leg = ['rms=' num2str(resnorm) ' ' leg] ;
        legen(cnt,1:length(leg))=leg;
    end
end
%figure ; h=plot(x,y) ; legend(h,legen);
alp = params(end,1) ; 
mug = params(end,2) ; 
lan = params(end,3) ; 
ome = params(end,4) ; 
fitted = y(:,end) ; 
end



function F = myfun(x,xdata)
F   = x(1)*(1-exp(-x(2)*xdata))+(1-x(1))*0.5* (1+erf((xdata-x(3))/(x(4)*sqrt(2))));
end






% alp= linspace((0),(1),5) ; % [0 1]
% mug= linspace((0),(1),5) ; % [0 1]
% lan= linspace((0),(1000),5) ; % >0
% ome= linspace((0.00001),(1000.00001),5) ; % >0
% 
% alp= 1-logspace(log10(0.00001),log10(0.9999999),5) ; % [0 1]
% mug= logspace(log10(0.000001),log10(0.9999999),5) ; % [0 1]
% lan= logspace(log10(0.01),log10(1000),10) ; % >0
% ome= logspace(log10(0.01),log10(1000),10) ; % >0
% 
% if exist('apriorimug','var')~=0 ; mug= max([0.01 apriorimug-0.2]):0.04:min([0.99 apriorimug+0.2]) ; end
% if exist('aprioriome','var')~=0 ; ome= max([0.01 aprioriome-aprioriome/5]):aprioriome/40:aprioriome+aprioriome/5 ; end
% if exist('apriorilan','var')~=0 ; lan= max([0.01 apriorilan-apriorilan/5]):apriorilan/20:apriorilan+apriorilan/5 ; end
% 
% cnt = 0 ; 
% pts = 0 ;
% rms = 0 ; 
% params = [0 0 0 0] ; 
% x = zeros(length(xdata),length(alp)*length(mug)*length(lan)*length(ome));
% x = zeros(length(xdata),1);
% y=x;
% chi2m=1000;
% lbo = [0;0;0;0.00001];
% ubo = [1;1000;1;1000.00001];
% options.Display='off';
% options.Diagnostics='off';
% options.FunValCheck='off';
% 
% for a=alp
%     for m=mug
%         for l=lan
%             [Ycalc,leg] = cumnomalizedcreationrate(Xobs,a,l,m,ome) ;
%             for o=ome                
%                 [par,resnorm] = lsqcurvefit(@myfun,[a l m o],xdata,ydata,lbo,ubo);
%                 
%                 if (resnorm<chi2m)
%                     chi2m=resnorm ;
%                     cnt=cnt+1;
%                     params(cnt,1:4) = par([1 3 2 4]) ; %[a m l o] ;
%                     
%                     x(1:length(xdata),cnt)=xdata; 
%                     [y(1:length(xdata),cnt),leg] = cumnomalizedcreationrate(xdata,par(1),par(2),par(3),par(4)) ;
%                     leg = ['rms=' num2str(resnorm) ' ' leg] ;
%                     legen(cnt,1:length(leg))=leg;
%                 end
%             end
%         end
%     end
% end
% whos x y legen
% figure ; h=plot(x,y) ; legend(h,legen);
% alp = params(end,1) ; 
% mug = params(end,2) ; 
% lan = params(end,3) ; 
% ome = params(end,4) ; 
% fitted = y(:,end) ; 
% end
% 
% 
% 
% function F = myfun(x,xdata)
% F   = x(1)*(1-exp(-x(2)*xdata))+(1-x(1))*0.5* (1+erf((xdata-x(3))/(x(4)*sqrt(2))));
% end





% % Call fminsearch with a random starting point.
% start_point = [0.5 0.5 5 5 ];
% 
% if exist('apriorimug','var') == 1 ;     start_point(2) = apriorimug ; end
% if exist('apriorilan','var') == 1 ;     start_point(3) = apriorilan ; end
% if exist('aprioriome','var') == 1 ;     start_point(4) = aprioriome ; end
% 
% model = @expfun;
% estimates = fminsearch(model, start_point);
% % expfun accepts curve parameters as inputs, and outputs sse,
% % the sum of squares error for A*exp(-lambda*xdata)-ydata,
% % and the FittedCurve. FMINSEARCH only needs sse, but we want
% % to plot the FittedCurve at the end.
%     function [sse, FittedCurve] = expfun(params)
%         A = params(1);
%         MU= params(2) ;
%         lambda = 10 ; params(3);
%         omega =  params(4);
%         % FittedCurve = A .* exp(-lambda * xdata);
%         FittedCurve = A*(1-exp(-1*lambda*xdata)) + (1-A)*0.5*(erf((xdata-MU)/(omega*sqrt(2)))) ;
%         ErrorVector = FittedCurve - ydata;
%         sse = sum(ErrorVector .^ 2);
%     end
% 
% [sse, FittedCurve] = model(estimates);
% alp = estimates(1);
% mug = estimates(2);
% lan = estimates(3);
% ome = estimates(4);



% alp= linspace((0.00001),(0.9999999),5) ; % [0 1]
% mug= linspace((0.00001),(0.9999999),5) ; % [0 1]
% lan= linspace((0.01),(1000),10) ; % >0
% ome= linspace((0.01),(1000),10) ; % >0
% 
% % alp= 1-logspace(log10(0.00001),log10(0.9999999),5) ; % [0 1]
% % mug= logspace(log10(0.000001),log10(0.9999999),5) ; % [0 1]
% % lan= logspace(log10(0.01),log10(1000),10) ; % >0
% % ome= logspace(log10(0.01),log10(1000),10) ; % >0
% 
% if exist('apriorimug','var')~=0 ; mug= max([0.01 apriorimug-0.2]):0.04:min([0.99 apriorimug+0.2]) ; end
% if exist('aprioriome','var')~=0 ; ome= max([0.01 aprioriome-aprioriome/5]):aprioriome/40:aprioriome+aprioriome/5 ; end
% if exist('apriorilan','var')~=0 ; lan= max([0.01 apriorilan-apriorilan/5]):apriorilan/20:apriorilan+apriorilan/5 ; end
% 
% cnt = 0 ; 
% pts = 0 ;
% rms = 0 ; 
% params = [0 0 0 0] ; 
% x = zeros(length(Xobs),length(alp)*length(mug)*length(lan)*length(ome));
% y=x;
% for a=alp
%     for m=mug
%         for l=lan
%             %[Ycalc,leg] = cumnomalizedcreationrate(Xobs,linspace(0,1,10),linspace(0,1000,10),linspace(0,1,10),ome) ;
%             
%             for o=ome                
%                 cnt=cnt+1;
%                 %disp(['test ' num2str(cnt) '/'  num2str(length(alp)*length(mug)*length(lan)*length(ome))])
%                 Ycalc = a*(1-exp(-l*Xobs)) + (1.0-a)*0.5*(1+erf((Xobs-m)/(o*sqrt(2)))) ;
%                 
%                 inds = logical(isnan(Ycalc)) ;
%                 Ycalc(inds) = -1;
%                 inds = logical(1-isnan(Ycalc)) ;
%                 if sum(inds) == length(inds)
%                     pts(cnt) =  sum(inds) ;
%                     rms(cnt) = sum((Yobs(inds)-Ycalc(inds)).^2)/pts(cnt) ;
%                     params(cnt,1:4) = [a m l o] ;
%                 else
%                     cnt=cnt-1;
%                 end
%                 leg=['alp=' num2str(a) ' mug=' num2str(m) ' lam=' num2str(l) ' ome=' num2str(o) ];
%                 %if nanmin(Ycalc) <0.1 ; disp(['nanmin(Ycalc) <0 for ' leg]);end
%                 
%                 x(1:length(Xobs),cnt)=Xobs; y(1:length(Xobs),cnt)=Ycalc; legen(cnt,1:length(leg))=leg;
% 
%             end
%         end
%     end
% end
% %figure ; h=plot(x,y) ; legend(h,legen);
% [RMS,ind] = min(rms);
% alp = params(ind,1) ; 
% mug = params(ind,2) ; 
% lan = params(ind,3) ; 
% ome = params(ind,4) ; 
% fitted = alp*(1-exp(-1*lan*Xobs)) + (1-alp)*0.5*(erf((Xobs-mug)/(ome*sqrt(2)))) ; 

